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PACS 47 . 20 . Ma - Interfacial instabilities (e.g., Rayleigh-Taylor) 
PACS 47 . 70 . -n - Reactive and radiative flows 
PACS 47 . 45 . Ab - Kinetic theory of gases 

Abstract -Reactive Rayleigh-Taylor systems are characterized by the competition between the 
growth of the instability and the rate of reaction between cold (heavy) and hot (light) phases. 
We present results from state-of-the-art numerical simulations performed at high resolution in 2d 
by means of a self consistent lattice Boltzmann method which evolves the coupled momentum 
and thermal equations and includes a reactive term. We tune the parameters affecting flame 
properties, in order to address the competition between turbulent mixing and reaction, ranging 
from slow to fast-reaction rates. We also study the mutual feedback between turbulence evolution 
driven by the Rayleigh-Taylor instability and front propagation against gravitational acceleration. 
We quantify both the enhancement of flame propagation due to turbulent mixing for the case of 
slow reaction-rate as well as the slowing down of turbulence growth for the fast reaction case, 
when the flame quickly burns the gravitationally unstable phase. An increase of intermittency at 
small scales for temperature characterizes the case of fast reaction, associated to the formation of 
sharp wrinkled fronts separating pure hot/cold fluids regions. 



7—< ' Introduction. — Many natural and industrial pro- 
^ oesses involve fluid transport and mixing of passive or 

• ^ dctive scalar fields; examples include concentration fields 

^^of chemicals or biological species as well as the temper- 
j^dture field in natural convection. While many of these 

■ ■ phenomena have been the subject of in-depth studies, the 
cases where chemical reactions are involved, presenting 
an even richer phenomenology, have received considerably 
less attention. We address the problem of the interplay 
of reaction and turbulent mixing in Rayleigh-Taylor (RT) 
systems (a situation occurring, for example, in thermonu- 
clear burning of type la supernovae [IHS] or in the iner- 
tially confined nuclear fusion [3]) focusing on the different 
regimes which develop as we vary the ratio between the 
characteric time scales of underlying turbulence, Tturb, and 
the reaction time, tr. We limit this study to the case of 



single-step reaction, i.e. the two reactant scalar fields are 
distinguished by a reaction progress variable, proportional 
to the temperature (see Fig. [1]). The two different temper- 
atures in the hot and cold blobs of fluid of our numerical 
setup, mimick the combustion of a cold mixture of ac- 
tual reactants into a hot mixture of burnt products [SHT] . 
The interesting point in this set-up is given by the natural 
competition between gravitational forces, which tends to 
mix the fluid and to produce a larger and larger mixing 
layer with uniform temperature, and combustion, which 
works against this mixing, trying to burn the whole vol- 
ume and producing a propagating flame of given tickness 
and velocity. Moreover, the global phenomenology is com- 
plicated by the natural unsteadiness of the underlying RT 
problem. The Damkohlcr number. Da, is the natural con- 
trol parameter and is identified by the ratio between the 
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turbulent time scale, Tturb and the reaction time scale tr- 
Notice that because of RT unsteadyness the Da number 
depends on time: 

Da{t) ^ Tturb{t)/TR. 

where Tturb oc t, as of standard RT phenomenology [7]. 

We perform highly resolved numerical simulations in 2d, 
with a resolution up to 4096 x 10000 grid points (see table 
[T]). The 2d set up allowed us to reach a wide scale separa- 
tion and a time-span large enough to address problems at 
both small and large Damkohler numbers, something still 
unfeseable in 3d. Our study has also a methodological 
motivation. We adopted a fully consistent thermal lattice 
Boltzmann method to evolve simultaneously the momen- 
tum equations and the advection-diffusion-reaction equa- 
tions for temperature. The novelty here is to show that 
the method works well also in a non trivial case where 
the thermal modes arc directly forced by the combustion 
terms. 

The main result of the paper concerns the quantification 
of the front propagation due to turbulence enhancement 
for the slow reaction case. Da <C 1, and the clear signature 
of a strong feedback on the fluid evolution induced by the 
flame propagation when Da > 1. In the latter case, we 
also measure an important increase of the temperature 
intcrmittcncy at small scales. 

Equations of motion and numerical setup. We 

adopt a numerical scheme based on a recently proposed 
thermal lattice Boltzmann algorithm [Hlin], which is able 
to reproduce the correct thermohydrodynamics of an ideal 
gas with good numerical accuracy [T^. To do that, the 
probability densities fi{x,t) for a particle with velocity 
ci (belonging to a discrete set, with the index I running 
over 37 values [3) at space location x and time t evolve 
according to the lattice Boltzmann BGK equation [TIHT5] 

Mx + ciAt,t + At) - fiix,t) = -—U- Z/'^')) {x,t)- 

TLB ^ ' 

(1) 

the Ihs stands for the free streaming of particles and 
the rhs represents the relaxation process towards Maxwell 
equilibrium //^''■'(x,^) with a characteristic time t^b (Ai 
is the simulation time step). Once the density (p), velocity 
(ti) and temperature (T) flelds are defined in terms of the 
lattice Boltzmann populations as 

p = Y^ fi; pu = Y, fici- DpT = Y,fi \ci - «l' , (2) 

(_D is the number of space dimensions), it has been shown 
[SlU] that the following set of macroscopic equations can 
be recovered (repeated indexes are summed upon): 

Dtp = -pdiUi 

pDtUi = -dip - pgSi^2 + vAui (3) 
pCuDtT = kAT + ^R{T), 
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Fig. 1: Initial configuration for the Rayleigh- Taylor system 
with combustion: cold fluid (fresh fuel) at T = on top and 
hot fluid (burnt material) at T = 1 on bottom. Such tempera- 
ture jump at the interface is smoothed by a hyperbolic tangent 
proflle with a width of the order of 10 grid points and with 
a randomly perturbed centre (thus enabling to perform inde- 
pendent runs). The system used is two dimensional and has 
size Lx X Lz plus periodic boundary conditions applied in the 
streamwise (x) direction. The fluid used is an ideal gas. 



where k is the thermal conductivity, tr the reaction time, 
g the gravity and Cy = D/2 the specific heat at constant 
volume. Eqs (jS]) are obtained provided that the thermo- 
hydrodynamic fields appearing in the equilibrium density 
functions are properly shifted: 

u^u + tlbO; T -fT + ^— '-g^ + R{T). 

D TR 

^ ^ (4) 

The novelty here is in the extra reaction term introduced 
for the temperature field 0. 

In the third equation of ([3]) we have already subtracted 
the compression term pV • u to avoid effects due to a vary- 
ing heat capacity or global heating of the system coming 
from steady increase of the underlying mean pressure. 

The reaction rate must be zero in the pure phases, which 
we set at temperatures T = 1, for the hot fluid at the bot- 
tom, and r = 0. for the cold fluid on top (see flgurc[T|), 
so that i?(0) — R{1) — 0; it must also transforms, irre- 
versibly, the pure cold phase (unstable) into the hot one 
(stable). A simple model for R{T) with these properties 
is given by a logistic-type expression 

R{T)^T{1-T) 

originally proposed [HlfTS] as a model for the propagation 
of an advantageous gene in a population. 

^This shift represents a kind of implicit equation, since, in prin- 
ciple, -R(r) should be a function of the "real" thermodynamic tem- 
perature, which must be shifted itself [3]; however, we can overcome 
this problem, observing that for th 2> t^b (always true in actual 
situations), and since the other shift is 0{g'^) (hence 0{Kv?) |9]) 
we can safely assume that R{T^^'>) fl(T). 
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Lx 




V 


9 




Tu 


Td 




r 


run A 


0.05 


4096 


10000 


0.005 


2.67 X 10 


-5 


0.95 


1.05 


5 X 10^ 


5.5 X 10* 


run B 


0.05 


4096 


10000 


0.005 


2.67 X 10 




0.95 


1.05 


5 X 10* 


5.5 X 10^ 


run C 


0.05 


4096 


10000 


0.005 


2.67 X 10 


-5 


0.95 


1.05 


5 X 10^ 


5.5 X 10^ 



Table 1: Parameters for the three sets of runs. Atwood number, At = {Td — Tu)/{Td + Tu); viscosity u (thermal diffusivities are 
the same since Prandtl number is 1 for each run; gravity g; temperature in the upper half region, T^; temperature in the lower 
half region, Td; reaction characteristic time tr; normalization time, r — \J Lx/ {g At). 



We performed three high resolution sets of runs (on lat- 
tices of 4096 X 10000 grid points) on the QPACE Super- 
computer |16lll7j . with different reaction times (run pa- 
rameters are collected in table [ij. For each set, we carried 
out several (0(10)) independent runs, in order to enhance 
statistics. 

Results and discussion. Any RT system, even in 
the case of r^j ^ 1, will eventually reach the fast reaction 
limit, i.e. a situation where Da{t) » 1. This is due to 
the fact that the underlying turbulence slows down adi- 
abatically, T{t) oc t. As a consequence, sooner or later 
the flame tends to become active, burning at a rate faster 
than the turbulence stirring/mixing. Here we study the 
two regimes Da ^ 1 and Da > 1 and the transition be- 
tween them. 

Mean temperature profiles evolution. For large Da, 
the mixing is effective only at very small scales (where the 
characteristic times of the fluid motion are shorter), while 
the reaction tends to make uniform the mixed regions: as 
a result we get a topology of the temperature field which 
is made of "patches" separated by rather thin interfaces, 
which are smoother than the non-reacting RT case [7|; 
in addition, the front of the hot phase moves, on average, 
with a non zero mean drift velocity towards the top. These 
preliminary features can be better understood, at a picto- 
rial level, looking at figure^ where we show the magnitude 
of the temperature gradient |Vr|^ = {{d^T)'^ + (d^T)'^) at 
three different times in the evolution for the fastest reac- 
tion rate that we have studied (top panel), and compare it 
with the non- reacting case (bottom panel). On the other 
hand, the larger the reaction time Tfj the closer is the 
phenomenology to the standard RT case: to see this we 
compare in figure[2]the evolution of the mean temperature 
profile 

f{z,t) = ^ I T{x,z;t)dx (5) 
Lx J 

for the two extreme cases in our database, runs A and C: 
while for rj^ = 5 x 10^ the evolution is basically undis- 
tinguishable from the usual RT dynamics [21 HH]: in the 
fast rection case (t/j = 5 x 10^) the center of mass of the 
system clearly moves upwards, due to the burning pro- 
cesses, causing a shift -and an asymmetry- of the mixing 
region. The propagation of the burnt hot material front 
against the fresh reactant (T = 0) can be quantified by 
the barycentric coordinate Zf{t), that we define as the 




Fig. 2: Snapshots of the magnitude field of the temperature 
gradient for the fast reaction case, run A (top panel), and for 
the non-reacting case (bottom panel). 



following integral [T51IT5] : 

Zf{t) ^ / T{z,t)dz (6) 

in figure m we plot the function Zf{t) vs t for the three 
different reaction rates: the growth of Zf{t) in time is 
greatly enhanced when going from small to large tr. 

Front propagation speed. If we integrate eqs. ([3]) over 
the whole volume, and divide by Lx, we get an exact equa- 
tion for the propagating front speed: 

Vf[t) = dt( r'^''\{z,t)dz\^-{T{l-T)), (7) 

(where ((•)) = {l/Lx)J J{-)dxdz) since the boundary 
terms vanish, owing to the periodic conditions on the lat- 
eral walls and to the adiabatic condition at top and bottom 
plates (VT|2^±2.j/2 = 0). For the laminar flame (that is 
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Fig. 3: Mean temperature profiles at various times for run A 
(bottom panel) and run C (top panel). The latter case is almost 
identical to the non-reacting case. 



without gravity producing turbulence), the integral can 
be evaluated exactly (using, for instance, the usual hyper- 
bolic tangent profile) to give an explicit expression for the 
speed, that is Vf oc Lf/rn, where i/ is the flame thick- 
ness: as the latter can be estimated to be L/ oc ^kt^, we 
end up with the well known result: 

Vf (X , (8) 

that is the flame propagates at constant speed. 
We now ask what changes when turbulence sets is. In 
the small Da limit, when turbulence has the time to mix 
the fluids before reaction becomes active, we are in the 
so-called pre-mixed combustion. In this case, it has been 
conjectured [THIIIO] that the simplest way to extend the 
result of the laminar case is to replace in expression (jS]) the 
molecular diffusivity k with an effective (turbulent) eddy 
diffusivity kt- If we use for the latter the dimensional 
estimate Kxit) ~ U{t)L{t), where U and L are large scale 
characteristic velocity and length (in our case, e.g. the 
root mean square velocity and the mixing region length), 
and plug it into ([8]), we get: 

V,it)^,[^^,^^^Uit),^^^lIWl^ (9) 

\ Tr \j Tr \j Tr 

where 

Vf{t) ^ U{t).l^ ^ U{t) Da{ty/\ (10) 



This prediction, probably valid to describe the evolution 
of slow flames in stationary turbulent flows is unlikely to 
be relevant for RT turbulence. The reason is that in or- 
der to observe a "eddy-diffusivity" driven flames one needs 
also scale separation between the turbulent eddies and the 
flame tickness, something that is not realized by the evolv- 
ing RT system. On the other hand, we can rewrite ([7]) 
exactly as: 

Vf{t)^-mi-f))^{W)]. (11) 

Tr 

where with 9 = T — T we denote the fluctuations with 
respect to the mean vertical profile. It is clear now that 
for Da < 1, the flame cannot have any strong influence 
on the underlying RT evolution and we can identify the 
first term on the rhs as the mixing layer length L{t) = 
(T(l — T)). Moreover, we know that in RT temperature 
fluctuations are almost constant in time and homogeneous 
inside the mixing layer [3T], so also the second term on 
the rhs is proportional to the mixing layer extension. A 
natural prediction for Da < 1 is therefore: 

Vf{t) cx — ; Vf{t) cx U{t)Da{t) (12) 
Tr 

In standard (stationary) turbulent reacting systems, one 
can check this prediction against experiments/simulations 
at various Da, obtained changing the reaction rate or the 
underlying turbulent intensity, while in our reacting RT 
setup we can exploit the fact that Da = Da{t) varies in 
time. In figure [S] we plot the front speed (normalized with 
the root mean square velocity) as a function of Da (which 
is itself a function of the simulation time) for the three 
runs. As one can see, our prediction (jl2p works satisfac- 
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Da 

Fig. 5: Front speed normalized by the root mean square (ver- 
tical) fluid velocity for the three runs as a function of the 
Damkohler number Da{t). The solid line represents the the- 
oretically predicted behaviour Vf/U oc Da, obtained on the 
basis that for Da < 1 flame propagates inside the well mixed 
mixing layer. The prediction Vf/U oc Da^^^, obtained from 
the assumption that in the pre-mixed combustion (slow reac- 
tion) regime one can simply substitute the molecular diffusivity 
with the turbulent one ft ^ kt in the expression for the lami- 
nar flame speed, is also plotted (dashed line). 



torily in a wide range of Da{t), showing deviations only 
for very small times, where turbulence is not yet devel- 
oped and the flame evolution is strongly influenced by the 
initial configuration, and for Da{t) > 1 where it cannot 
be expected to be valid. In the latter case our data point 
flatten, as we clearly observe the feedback of the flame 
on the turbulent evolution, with a sort of synchronization 
between flame propagation and evolution of the turbulent 
kinetic energy toward a value where Vf{t) U(t). Such a 
behaviour turns out to be in agreement with recent the- 
oretical results obtained through a mean-field approach 

Small scale intermittency. When the reaction rate is 
fast {Da ^ 1), there are no extended regions which are 
well mixed, since the cold material, as soon as it is slightly 
entrained through the hot one is suddenly burnt. As a re- 
sult, the temperature field organizes in patches of pure 
reactants/products separated by sharp interfaces (being 
in the so called "segregated regime"), and, consequently, 
it has been conjectured that an increased intermittency 
develops at the small scales [7]- The authors in [7] also 
derived a phenomenological prediction for the scaling laws 
of fluid temperature (and velocity) structure functions, ac- 
cording to which, in the asymptotics of Da ^ 1, they 
should follow the relation 

S^^\R,t)^{\SHm^(^-^y\ (13) 



1 I , , 1 

10' Iff- 10" 10"' 

m/ia At) 

Fig. 6: The 4-th order flatness F4, for the three runs and for the 
non-reacting RT {tr — > 00). Data from run A {Da 3> 1) agree 
well, within error bars, with the prediction given by equation 



(where L{t) is the mixing length), irrespective of the order 
p. From Eqn. (|13p the expression for the flatness reads: 

and so it increases with decreasing R for all orders, a clear 
indication of strong small scales intermittency. In figure IHl 
we show the growth of ^4 as function of the mixing length 
L, for the three runs: the flatness for run A, corresponding 
to the smallest reaction time, is in good agreement, within 
error bars, with the prediction of equation Ft ~ 

2^2/3. ijistead, at increasing Tij, intermittency is depleted 
and the flatness grows more slowly, at a rate comparable 
(within error bars) with the non-reacting RT case, whose 
data are also reported for comparison. 

Conclusions. — We used a self consistent thermal 
lattice Boltzmann algorithm to perform numerical sim- 
ulations of 2D Rayleigh-Taylor turbulence, in presence of 
chemical reactions between hot and cold fluids. The re- 
action was modelled by means of a Fisher-Kolmogorov- 
Petrovsky-Piskunov source term in the temperature equa- 
tion; this term has been introduced by application of a 
suitable shift of the temperature field appearing in the 
equilibria of the lattice Boltzmann equation. 

We analyzed the crossover among the various regimes 
emerging from the competition of turbulent mixing and 
reaction, going from the segregated {tr Tturt) to the 
well mixed one {tr ^ Tturb)- We showed that, in the 
latter case, the effect of turbulence is to enhance the re- 
action front speed leading to an homogeneous burning in 
the whole mixing layer region. On the other hand, for 
moderate and large Damkohler, there is a feedback of the 
reaction on the statistical properties of the temperature 
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field, resulting in increased intermittency at small scales 
in reasonable accordance with the prediction of [7] . 
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